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ABSTRACT 

Since the original work of Hansen and Tisserand in the XlXth century, there have been many variations in the analytical 
expansion of the three-body disturbing function in series of the semi-major axis ratio. With the increasing number of 
planetary systems of large eccentricity, these expansions are even more interesting as they allow us to obtain for the 
secular systems finite expressions that are valid for all eccentricities and inclinations. We revisited the derivation of the 
disturbing function in Legendre polynomial, with a special focus on the secular system. We provide here expressions of 
the disturbing function for the planar and spatial case at any order with respect to the ratio of the semi-major axes. 
Moreover, for orders in the ratio of semi-major axis up to ten in the planar case and five in the spatial case, we provide 
explicit expansions of the secular system, and simple algorithms with minimal computation to extend this to higher 
order, as well as the algorithms for the computation of non secular terms. 

Key words. Celestial mechanics-N-body problems-Planetary systems-Methods: analytical 



1. Introduction 

In the three-body problem, there are two classical ways 
to compute the principal part of the disturbing function 
—Qmm'/A. The first approach is to expand it in powers of 
eccentricity and inclination, with coefficient s that are ex- 



pressed in term of Lapla c e coefficients (e.g. Laplace 1785t 



I Abu-El- Ata fc Chanrontl 119751: lLaskar fc Robutell 119951 
but this approach, which is well suited to the study of 
the Solar System, has some limitations for some extra-solar 
planetary systems, where the eccentricity can reach very 
high values. Another drawback, is that without some spe- 
cial truncation corrections, the angular momentum will not 
be conserved exactly in the truncated system. 

In the second approach, the expansion is made with re- 
spect to the ratio of the semi-major axes of the two bodies 
a = a/a' , where a' is related to the external body. Although 
it may be less efficient for low eccentricity and large val- 
ues of a, as for the inner planets in the Solar system, the 
advantage is that this expansion allows us to obtain fi- 
nite expressions for arbitrary eccentricities in the secular 
system, while both approaches allow expansions for arbi- 
trary inclinations. The most important contr ibution to this 
problem was made in th e XlXth century ([Hansen! 118551: 
iHiill [18751: Ffisserandi n"899) . With the development of space 
technology, expansions of the disturbing function in term of 
Legendre polynomials have been rejuvenated for the con- 
struction of satellite theories in the vicinity of the Earth 
(lKozailll959t lKaulalll962t lBrumberg||1967t iBrumberg et al.l 
119711: lGiacaglialll974UAbu-El- Ata fc Chaprontlll975lf ~ 



The discovery of new plane tary systems, starting with 
51Pegb ()Mavor fc Quelodfl"9 95). has rasemi-majorised the 
need to revise these methods as many systems have plan- 
ets with very high eccentricity, as GL581, HD217107, 
HD69830, HD74156, HD168443, HD102272, HD169830, 
HD202206, HD183263, or even HD80606, where the ec- 



centricity reaches 0.933- With the discovery of these nu- 
merous new planetary systems, the previous analytical ex- 
pansions in Laplace coefficients developed for the Solar 
System are no longer the most appropriate, and we are 
faced with the need to understand more globally the dy- 
namics of these systems. On the other hand, as the pa- 
rameters of these extrasolar planetary systems remain not 
very well known, there is no necessity for very precise ana- 
lytical approximations. This has led several authors to use 
the Legendre expansion of the potential for the dynam- 
ical study of the secular planetar y system, following the 
previous studies of stellar systems ( Krvmolowski fc Mazehl 
119991: iFord et al.l l2000t iBlaes et all 120021) . where the secu- 
lar spatial three-body system wa s expanded up to th e oc- 
tupolar order (a 3 ). In particular. iLee fc Peald ([20031 ) used 
the secular planar system at octupolar order to study the 
secular dynamics of the HD 1 68443 and HD12661 systems. 
Migaszcwski & Gozdzicwski (2008) computed the planar 
secular system to high order using computer algebra to av- 
erage over the mean anomaly. 

As it appears that there is a growing interest in analyt- 
ical studies of extra solar planetary systems, we considered 
it interesting to present a derivation of the planetary (or 
stellar) disturbing function in a very simple and explicit 
way that does not already appear in the existing literature. 
We aimed to write a self-contained paper that allows one 
to construct explicitly the planetary Hamiltonian to high 
order, with minimal additional computation. In particular, 
in Sect. [3] we show that the planar secular system can be 
obtained explicitly at any order without the need for com- 
puter algebra. The spatial case is treated in Sect. @]with 
minimal computation when expressed with respect to the 
mutual inclination J. In the presence of more than two 
planets, it may be more convenient to use a fixed reference 



see http://exoplanet.eu 
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frame, and the derivation of the spatial Hamiltonian is pro- 
vided in this case in Sect. 14.21 Our deriva tion is close to 
the or iginal formulations of H ansen! (|1855| ) and iTisseraiidl 
(1899), but with emphasis on the secular system and the 
direct derivation of explicit expressions. We also present in 
the Appendix [A] a new derivation of the Hansen coefficients 
for the secular terms. 



2. Expansion of the disturbing function 

To simplify the notations, and although everything can 
be generalized to an TV-body system, we consider here a 
three-body problem with a central body of mass M and 
two other celestial bodies of masses m, to'. If we write the 
Hamiltonian of the Newtonian interactions among these 
three bodi es in Poincare canonica l heliocentric variables, 
we obtain (jLaskar fc Robutellll995l Eq. 15) 



where 



'Ho — — 



G(M + to) G(M + to') 



2 a 



2a' 



(2.1) 



(2.2) 



is the Keplerian interaction with elliptical elements 
a, e, i, M, uj, and f2 that denote, respectively, the semi-major 
axis, eccentricity, inclination, mean anomaly, argument of 
perihelion, and longitude of the node (with primes for the 
external body of mass to'). The principal part IA\ part of 
the perturbation and indirect part T\ are then 



3. Planar case 

We first study the planar case that leads to some simpli- 
fications in the expansions. We define v, v' to be the true 
anomalies of u, u', and uj,uj' their argument of perihelion, 
u = v + uj, u' = v' + uj', and x = u — u' . We have 
then u-u' = cos a: and, followin g a classical computation 
fe.g. IWhittaker fc Watsor] fl927l p. 303), we have for all 

*e[o,i[, 



(l-2zcos2;-r-z 2 )- 1/2 = (1- ze lx )-' r \l- ze- lx ) 



1/2, 



EE 

p=0 q=0 



(2p)\{2q)\ 



22p+2g(p!)2(g!)2 



e l(p-q)x z p+q 



and, with n = p + q, and after changing q to n — q, 
{2q)\{2n-2q)\ (2q _ n)x \ 



Thus 



with 



Sv§ 22n (« ! ) 2 «"-9)!) 2 * 



± n / . J n,q^ 
q=0 



fit 



(2q)\(2n-2q)l 



-1/2 

(3.1) 
(3.2) 

(3.3) 
(3.4) 



' n ' q 2 2 ™(g!) 2 ((n-g)!) 2 

for < q < n. If we write x — v ~ v' +lo — lu' , T n can now 
be expressed in the form 



-G- 



Ti 



M 



(2.3) 



where r, r' are the radius vectors of the inner and outer 
planets, with norms r,r', unit vectors u = r/r and u' = 
r'/r', and conjugate momentum f,f. We focus first on the 
principal part of the Hamiltonian Hi, which is the most 
difficult part to compute, while the computation of the in- 
direct part will be made in Sect. [5J With p = r/r' , and 
T = a! I ||r — r'||, we have 



T= -(l+p 2 -2pu-u' 



-1/2 



/ oo 



-^P n (u-u'K, (2.4) 
n=0 



where P n {x) are the Legendre polynomials that can be writ 
ten as 

[n/2] 

p n {z) = y j>„., i n 



with 



Pn,k 



(2n - 2k)\ 



2™ k\{n - k)\{n - 2k)\ 



(2.6) 



With a = a/a',^ — r/a."/' — r' /a', we have p = 07/7', 
and thus, as 7, 7', u, u' do not depend on a, a', the expan- 
sion of T in powers of a is 



o 



and 



with T n 



72] 



7 



,/n+l 



Fn 



F n = P n (u-u') = J2Pn,k(n-u') n - 



2 A- 



(2.7) 



(2. 



k=0 



J~ n ^ J~ n.q & 

q=0 



(2q—n)(uj—uJ ) 



with 



J^n.q — Jn,q m+l 



2. e 1 ( 2 9-")( t '- 1 '') 



(3.5) 



(3.6) 



The quantities F n q can then be expressed in term of 
Hansen coefficients X^ ,m (e) defined for n, to G Z as 



+00 

(£j e imv = X'k' m (e)e lkM 



(3.7) 



k= — 00 



For convenience, we denote X T ^' m = X^' m (e) and X'2 ,m 
X^ m (e'). We have thus 



+ OO 



(2-5) F n , q = f n , q £ X l 



2q-n ^-/-(n+l),-2g+n ^(kM+k' 'M') 



k,k' = — 00 



(3.8) 

For arbitrary values of fc € Z, the Hansen coefficients X^' m 
can be expressed in an explicit manner as an infinite se- 
ries i nvolving Besse l functions and hypergeometric func- 
tions (|Hansenl 118551: lTisserand[ 18991). or in terms of gen- 
eralized Laplace coefficients ~(lLaskarll2005h . but for k = 0, 
X™' m reduces to a finite polynomial expression in e, 1/e, 

y/1 — e 2 , and 1/y/l — e 2 (see Appendix A). We have thus a 
very compact expression for the coefficient of any argument 
kM + k'M' in explicit form at all orders n in a. Indeed, if 

( k k' ) 

we denote this coefficient T n ' , that is 

+00 

F« = E J^ fc V( feM+fe ' M ') , (3.9) 
k,k'=— 00 
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F 


= 1 


Fi 


= cos a; 


F-2 


= (l + 3cos2a;)/4 


F 3 


= (3 cos x + 5 cos 3x) /8 


Fi 


= (9 + 20 cos 2x + 35 cos Ax) /64 


F 5 


= (30 cos x + 35 cos 3a; + 63 cos 5a; ) /128 


F e 


= (50 + 105 cos 2x + 126 cos 4a- + 231 cos 6a-) /512 


P- 

rj 


— i i i o lu» a ^ loy cos oa zoi lub ua ^ izy coo / a j / l\j£^± 


F 8 


= (1225 + 2520 cos 2a; + 2772 cos 4a; + 3432 cos 6a; + 6435 cos 8a;) /16384 


F 9 


= (4410 cos x + 4620 cos 3a; + 5148 cos 5a; + 6435 cos 7a; + 12155 cos 9a;) /32768 


F w 


= (7938 + 16170 cos 2a; + 17160 cos 4a; + 19305 cos 6a; + 24310 cos 8a; + 46189 cos 10a- ) /131072 



Table 1. Tisserand functions for the planar case (Ea l2.8[) . We have x — u — u' with u = v + to, u' = v' + ui' . 



we have which allows an immediate computation of any argument 

n 

-rVk.k') _ f y n,25-n y r(n+l),-2?+n i(2ii-ii)(uV) 



J-~n k ' k ^ in terms of Hansen coefficients. 



9=0 



-(0,0) 



(3.10) 3.2. Computation of the secular part J 7 ^ 

In particular, for all neN, the secular part J>r 0) can The computation of the secular part of the Hamiltonian 

6V ^ be Simplified ' using the classical relation X n j- m = jrjPfi) ig quite gimplej ag we wiU jugt haye tQ make the 

X^ ,m among Hansen coefficients and the relation f n>q — transformation 



/„,„_, from Eq.(03) y ,-(n+i),m , , 

cos(ma;) — > X ' X cos(m(cj — cj )) (3.14) 



in the expression of F n . Moreover, we have x'- (n+1) ' m = 



Tt(O.O) c f yii,0 7 /-(»+1),I), 
J-n — e nJn,%X X Q + 

[(n-i)/2] for m > n > 1 (see Appendix |A"| . The term in cosna; can 

V 2 f n , q X^ n ~ 2q X'o (" +1 )^- 2 « cos ((n - 2g)(w - u/)) > thus be discarded in the expression of F n (Table Q} which 

g=0 simplifies the expression of the Hamiltonian. The secular 

(3.11) Hamiltonian is thus expressed in finite form, using the val- 

where e„ = if n is odd, and e„ = 1 if n is even. ues of the Hansen coefficients given in Appendix E] We 

have 

3.1. Practical algorithm j-(o.o) _ 1 

Equations (|3.4[) and (13.111) provide an explicit algorithm J^ ' -* = 

for the computation of the principal part of the disturbing ( ,o) 1 2 o /-3.o 

function J 7 , and in particular for the computation of the ^"2 — 4^0 ^ 

secular Hamiltonian at any order N in a 

J _m 01 o 

iv 

<^>=^j-r ) «-, (3.12 

ra=0 



t-(0,0) v 3,l v /-4,l / /\ 

.F3 ' = -AT a o cos(w — cj ) 



-,-(0,0) 9 v 4,0t A /-5,0 5 v 4,2 v /-5,2 /\ \ 



for which the involved Hansen coefficients 

Xq ' m , X'q (" +1 ) ,m reduce to finite expressions (see J 7 ^ ' ^ = — Xq' 1 X'q 6,1 cos(w — lj') 

Appendix A). As there are no se cular terms in the indirec t ^4 

part of the Hamiltonian (see lLaskar fe Robutel 119951) . — Jf 5 ' 3 X'~ 6 ' 3 cos(3(a; lj')) 

the computation of the secular Hamiltonian reduces to 128 

the computation of the rational constants f 7h k given by , Q % 25 6 ,-70 105 6 2 ,-72 , 

Eq. p.4p . In practice, for finite order n, it is even easier to -^6 = 256^° ^ 512^° cos (2( w — w )) 
use the expression of F n p.3p (Table [T|) and to translate it 

into J 7 ^ (|3.11[) by the simple transformation _| Xq A X'q 7,4 cos(4(w — a/)) 



256" 

7 " T cos(mx) = jfo) = ^^'^'o" 8 ' 1 cos(. - c') 



7 ,n+i ' ' " ' 1024 

E v n.rn v ,-(n+l)-rn futfiU'Uf't I '\\ + TfYM^ ^ X ' ' C0S ( 3 ( W ^ W ')) 

X k X k , cos{kM + k M + m(uj — uj )) , 1U24 

k,k'=— oo 231 7 5 „/— 8,5 /_/ / \ \ 

(3.13) + 1024 X ° X ° cos ( 5 ( w -^)) 
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F = 1 

Fi = fi cos(m — v!) + v cos(u + it') 

1 3 3 3 3 3 3 

F2 = — g + 4^ + 4 A* 2 + ^M 2 cos(2u — 2m') + -^y 2 cos(2m + 2m') + -^uficos(2u) + ^u/icos(2u) 

( 3 15 2 15 3 \ / 3 15 3 15 2 \ 

F 3 =1 -^M + — f M + -g-A* I cos(m — u J + I — + -g-w + -^v[i j cos(m + u ) 

+ ^-v 2 [i cos(3m + u) + -^-u 2 'a cos(m + 3m') 
o 

15 15 5 5 

+ -^-ufi 2 cos(3m — u) + -5-vfi 2 cos(m — 3m') + „ u 3 cos(3m — 3m') + -^u 3 cos(3m + 3m') 

3 15 o 105 4 15 2 105 22 105 4 
=+ 8 ~ T ~64" ~ T M "16~ M ~64~ M 



15 2 105 
~8~ M + T6" 



"V + §M 4 ) cos(2m - 2m') + (~v 2 + + cos(2m + 2u') 



i 15 105 3 105 3 \ ln , / 15 105 3 105 , , 

+ ( — r J/ / i + v+ Te 11 ^ ) cos ( 2u ) + \ — 4-^ + -je"" m + "16"^ Jcos(2m) 
35 35 35 35 

+ cos(4m — 4m') + ttt^ 4 cos(4m + 4m') + -rr.vu? cos(4m — 2v!) + -r-vn 3 cos(2m — 4m') 
64 64 16 16 

105 99 / , \ 105 99 / . *\ 35 q , . _ /> 35 q *, 
H — gg" 1 ' A* cos(4m) + -^-v fi cos (4m ) + j^v /xcos(4u + 2w ) + ">cos(2m + 4« ) 

105 2 945 4 105 3 945 2 3 315 
-~8T U »+-6A v »-^6» +64' 

105 3 315 5 105 2 945 3 2 945 M 
-^g-f +"64"^ ~ ~8~^ + "32" !/At + ~6A 1/ ^ ) cos ( u + u ) 

2 315 4 945 2 3^ / n / \ ( 105 2 945 3 2 315 , > 
^ + "32"^ ^ + "64~^ ^ / cos ' u + 3u ) + I — ^[g"^/ 1 + "54"^ ^ + ~TI V ^ I cos ' u ~ 3u ) 
2 315 4 945 2 3 \ ,„ / \ I 105 2 945 3 2 315 A 
v u + -^2~v A 4 + "54"^ ^ / cos y^ u + K ) + ( ~~\q v L jl "64"^ ^ ~32 UfJ ' I cos y* u — u ) 

315 5 315 3 2 \ , n „ ,, / 35 3 315 2 3 315 

315 39 / K /\ 315 23 / k /\ 315 3 2 . 315 2 3 /\ 315 4 , /. 

H — ttt-v a cosim + ou ) + -rrrV u cosnt — 011 ) + -rrr-v a cos 5it + u ) + -rrrv U cos(5u — u ) + -r^v ucos(3u + 5m I 
64 64 64 64 128 

+ T^^u 4 cos(3m — 5m') + ^^^ 4 ucos(5m + 3m') + ^n^v^ 4 cos(5m — 3m') + S^v' J cos(5m + 5m') + r^-u 5 cos(5m — 5m') 
Izo 12o 12o 12o L2o 

Table 2. Tisserand functions for the spatial case. We have u = v + uj, u' = v' + ui' while fi — cos 2 ( J/2), v — sin 2 (J/2), 
where J is the mutual inclination. 



=1 


' 15 


+ 


( 15 


+ 


( 105 


"IT' 


+ 


( 105 


(-IB" 1 


+ 


( 35 

-is" 



These results are equivalent t o the expression obtained by 
Mig aszewski fc Gozdziewsk i (2008) using computer alge- 
bra. 

4. Spatial case 

The spatial case is more complicated as it involves addi- 
tional variables. Our goal is to derive explicit formulae that 
are as compact as possible. We thus expand in terms of 
the mutual inclination J. For each orbit, we use a reference 
frame (i,j,k) associated with the orbit, with first vector i 
in the direction of the ascending node of r' over r. With 
u = v + uj and vl = v' + a/, we have 

u-u' = coswcosm' + cos J sin u sin u' j\ 
= fi cos x + v cos y 

with the same notations as iTisseran 

x = U — v! ; v = u + u ; a = cos 2 — ; v = sin 2 — . 

' y ' p 2 2 

(4.2) 



As in the planar case, we have F n = P„(u-u'), but now u-u' 
is given by the slightly more complex expression (|4.1I) . For 
all n, we have 

n n 

F n =^2Yl Q s n v 0*> v) e lu («- 2s ) e " 1 ' (™- 2 «) , (4.3) 

s=0 g=0 

where the q!"^ (fJ>, v) are polynomials in fi, v of degree n that 
are call ed the Tiss e rand polynomial^ as a recognition of the 
work of lTisserandl (118851). althou gh these expressions are al- 
ready present in ((Hansen] [l 855). As (u,u') — > (—u,—u') 

leaves u-u' unchanged, we have Q^Z sn _ q — Qi?q, thus 
F n can be expressed as a trigonometric polynomial in 
cos(mw + m'u'). Although they can be computed explic- 
itly for all n (see Appendix B), for a given n, it is often 
more efficient to make a direct computation of F n on a com- 
puter algebra system. For example, -F20 is computed in less 

2 One should consult ( Aksenovl Il986t ) for a detailed discus- 
sion of the relation bet ween the Tiss erand polynomials and the 
inclination functions of iKaulal (|l962f ). 
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than 1.5 seconds in exact rational arithmetics with TRIP 
(|Gastineau fe Laskarl l2009) on an average laptop computer 
using the simple expression (12.81) . We can then express 
T n = (j n /j' n+1 )Fn in term of Hansen coefficients. We thus 
have 



+00 +00 



(kM+k'M') 



(4.4) 



k=—oo k'=—oo 



(0,0) _ v 4,0 w-5,0 



A o A 



105 , 15 , 

v u 2 

64 8 F 



105 , , 105 , 

v 2 n 2 H u 4 

16 P 64 P 



3 15 , 
1/ 

8 8 

^4,2 v /-5,2 
~ A A 

15 o 105 ri 35 j . , 

( u 2 H i/V H if ) cos(2w - 2u/) 

1 8 M 16 p r' v ; 



with 



, 15 o 35 a 105 00.. 
+ — ^ + ~" + -r^v V) cos(2w + 2w' 
8 16 16 



T {k,k') _ 

n 

En(") v™- n ~ 2s ■vW _ ( n + 1 )> n_2 9 ("-2s)w 1 (n-2«)w' 
^s,g A fc A fc' e e 

s,g=0 

(4.5) 

More practically, starting from the finite expressions in co- 
sine polynomials (Table [5]) of the Tisserand functions F n , 
the expression of J- n is obtained as in the planar case by 
means of the more general transformation 



15 

T 



105 
~16 



105 

+(~vfi + ^t^H + —vii a )X^X' u,u cos(2w) 



r4,2 v t-5,0 



15 105 r> 105 o x /in 1 — fi 2 , 

+(- T ^ + -jg-i^A* + — ^ 3 )a- 4 '°a" 5 ' 2 cos(2a/) 

. v 4,4 v /-5,2 

35 35 

H ^/x 3 cos(4o; - 2u/) H j/ 3 ^cos(4w + 2w') 

16 16 



105 
li2~ 



4,4 v/-5, 2,2 



X^X 







v n cos(4cj) 



-co +00 



j cos(mu + m'u') = 



(4.6) 



k— — ock'— — oo 

X%< m X'- {n+1) ' m ' cos(fcM + k'M' + mw + ro'w') . 



4.1. Computation of the secular part 

As in the planar case, the computation of the secular part of 

the Hamiltonian Fn is just a straightforward translation 
of F n in Table [2 with the transformation 



cos(mw + m'u' 



Moreover, as X 



' X cos(mw + m w ) 



/-(n+l), 





(4.7) 

for n > 1, all terms in 



cos(mu ± n«') can be discarded in F n . We thus have, 



77(0,0) _ 1 
7 -(o,o) = 



7". 



(0,0) 

1 3 



(-- + ^ + iM 2 )X 2 '°X' 3 '° + -^o"AV' U cos(2 W ) 



-2,2 v /-3,0 



T-(O-O) _ 

v 3.1 v /-4,l 
A A 



/ 3 15 9 15 , /x 

^2^ + T ^ + ~8~ M ^ cos ^ ~~ w ^ 



.3 15 o 15 n. , 

- (— -v + —v + —vu, ) cos(w + bJ ) 
2 8 4 

15 ,,3,3 w-4,1 r.,2 



+— X*' 3 X' • [u*n cos(3w + J) + !//z 2 cos(3w - J)) 



We thus observe here an important simplification in the 

quadrupolar secular Hamiltonian F^ ,0 \ as all terms in- 
volving the external planet longitude of perihelion oj' van- 
ish a nd we are left w i th an integrablc Hamiltonian. This is 



what Lidov fc Ziglinl ljl976l| called a happy coincidence (see 
also iFarago fc Laskarl 120101. This is no longer the case at 
higher orders. 



4.2. Expression in a fixed reference frame 

The above expressions are given with respect to the mu- 
tual inclination to shorten the algebraic expansions. This 
is especially useful in a three-body problem, but to obtain 
expressions valid in a fixed reference frame, then one needs 
to substitute into u-u' its expression in terms of the ellip- 
tical elements of the two bodies. We ca n then generalize 
the expression (14.11) and write it now as (Brumberg|[T967 
lAbu-El-Ata fc Chaprontlll975tlLaskar fc Robute]||199V 



u-u' = ^{^e lx + ^e ly ) (4.8) 
with as before x = u — v! , y = u + vf , and 



//* = {cc e 2 + ss e 
= ^cs'e 1 



(4.9) 



fi-sV . _ n — n / 

2 — see 2 



with c = cos(i/2), s = sin(i/2), and the same for the primes. 
We write 

/it* = a + ia' ; v* = b + ib' (4-10) 

with 



a = (c 2 c' 2 + s 2 s' 2 ) cos(n - f2') + 2cc'ss' 

a' = (cV 2 -s 2 s' 2 )sin(ft-ft') 

b = (c 2 s' 2 + s 2 c' 2 ) cos(n - n') - 2cc'ss' 

b' = (c 2 s' 2 -s 2 c' 2 )sin(^-fi') . 



(4.11) 
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F = 1 

Fi = —b' sin(it + u') + bcos(u + u) — a' sm(u — u) + a cos(u — u') 

f 2 =-i+^' 2 +^ 2 +r' 2 +!« 2 

33 33 33 33 

+ (^ba - -b' a') cos(2«) - (-6a' + ^b'a) sin(2u) + ( ~6'a' + ^ba) cos(2m') + (-ba' - ~ 6'a) sin(2u') 

+ (|& 2 - |fe' 2 ) cos(2u + 2m') + (|a 2 - |a' 2 )cos(2w - 2u') - |fc'fc sin(2u + 2u) + |a'a sin(-2u + 2u) 

,15 T /o 15,9 15*9 15 q 3, / / x .15./9 15,9 3, 15 , /o , 15 , ^ N , * x 

= +(—b ,2 a + — Jra + — a /2 a + — a" 3 - -a) cos(~u + u') + f-pfeo + -rha - -^b + -5-6 & + -5-6 ) cos(u + u') 
4 4 o o 2 4 4 z o o 

+ A' _ *jtys _ i^ 6 ' 6 2 _ IjV a ' 2 - ^fe'a 2 ) sin(u + u') + (^a'a 2 - jja' + ^&'V + ^-b 2 a' + ^a' 3 ) sin(-w + u') 
2 o o 4 4 o 2 4 4 o 

+ (-^&'6a' - ijVa + -^6 2 a) cos(u + 3u') + (^&V - ^-b' 2 a' - ^-b'ba) sin(u + 3u') 

4 a o o o 4 

+ (^fra 2 - ^6a /2 + —tfa'a) cos(-w + 3it') + (-5-6' a' 2 + -^-ba!a - -JVa 2 ) s\n(-u + 3u') 

o o 4 o 4 o ' 

o 4 o o o 4 

+ (^-ba 2 - ^6a /2 - -p&Va) cos(-3u + u') + (yWa + -^-b'a 2 - -5-&V 2 ) sin(-3u + u) 
o o 4 4 a o 

+(5& 3 - ^-b' 2 b) cos(3n + 3u') + (|fe' 3 - -^&'& 2 ) sin(3u + 3u') 

5 15 15 5 

+ ( o" 3 - -5" a ' 2a ) cos(-3u + 3u') + i^-a'a 2 - 3 a' 3 ) sin(-3u + 3m') 
o o o o 

Table 3. Tisserand functions for the spatial case in a fixed reference frame. We have u = v + u>,u' = v' + to' while 
a, a', 6, 6' depend on the nodes and inclinations and are given in (|4.11l) . 



With these notations, we have 

,2/ 



and 



and 



H = sin 2 (J/2), = cos 2 (J/2) 



u u' = a cos(u — v!) — a' sin(u — v!) 
+ b cos(u + v!) — b' sin(u + u') 



(4.12) 



(4.13) 



The remaining part is identical to the previous case. The 
expressions of the Tisserand functions F n are given for < 
n < 3 in Table [3] and for all n in Appendix [Bj In practice, 
for small values of n we can use explicit expressions of F n , 
using the straightforward approach consisting of computing 
F n directly with computer algebra, using the relation (|2.8[) , 
and then to translate it as previously in order to obtain the 
expression of any argument kM + k'M' by means of the 
relations 



7 



cos 

yn+1 \ gin 



(mu + m'u ! ) = 



k— — oc k'— — cx) 

x n, mxl -(n+l),m' j cos (fcM + # M> + ^ + m ^, } _ 

(4.14) 



We thus have for the secular part J 7 , 



(0,0) 



(0,0) 
3 

(0,0) 



F^ = X^X'- 3 '°[-i + \{b" + b- + a- + a 2 )] 



-3,0 



■[(6a - b'a') cos(2w) - (6a' + 6'a) sin(2w)] 



S3 - y 

4 

+(26' 2 a + 26 2 a + a' 2 a + a 3 - -a) cos(w - J) 

D 

4 

+(26a' 2 + 26a 2 - -6 + 6' 2 6 + 6 3 ) cos(w + uj') 
5 

+(-6' - 6' 3 - 6'6 2 - 26'a' 2 - 26'a 2 ) sin(w + a/) 
5 

4 

-(a'a 2 - -a' + 26' 2 a' + 26 2 a' + a' 3 ) sin(w - uj')] 
5 

(6 2 a - 26'6a' - 6' 2 a) cos(3a; + uj') 
+(6a 2 - 6a' 2 - 26'a'a) cos(3w - uj') 
+(6' 2 a' - 6 2 a' - 26'6a) sin(3w + uj') 
-(26a'a + 6'a 2 - 6'a' 2 ) sin(3w - a;')] ■ 

5. Indirect part 

Until now, we considered only the principal part of the per- 
turbing Hamiltonian. The computation of the indirect part 
7i (|2.3I) is more straightforward, as it compares with the 
computation of r ■ r' . Indeed, we have (jLaskar fc; Robutel 



15 



3,1 w—Mr 
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1199.1 Eq. 24), 



_ mm 

Ti = —-;-V with V = rv' 

1V1 



(5.1) 



where r, r' are the velocities in the corresponding Keplerian 
problem. With classical computation, we obtain 



a 2 

-cosi? = -> Jk(ke) cos kM 
r e ^ 

fc=i 

- sin E = 2 J 'k ( fce ) sin kM 

k=l 



(5.2) 



where E is the eccentric anomaly and Jk(x) are the Bessel 
functions. The coordinates (X, Y) of the velocity r in the 
reference frame of the orbit with origin at perihelion are 
then easily expressed in Fourier series of the mean anomaly, 
as 

'a 



X = 



-na | — s'mE 
r 



Y = na\/ 1 — e 2 [ — cosi? 



(5.3) 



If we denote Z = X + iY, we then have for the spatial prob- 
lem in the fixed referenc e frame, with the same n otations 
as in the previous section (|Laskar fc Robuteilll99l Eq. 37), 



i(u+a 



(5.4) 



If one considers the mutual inclination J, as in section |4j 
this expression simplifies since i7 = 57', s' = 0, c' = 1, and 
/x*,^* are then real, with //* = cos 2 (J/2), v„ = sin 2 (J/2). 
In the planar case (Sect. [3]), /i* — 1, — 0. 

We considered here heliocentric coordinates. One could 
also use Jacobi coordinates for the three-body problem. In 
this case, the indirect part does not require a dditional com - 
putations as it is expressed in term of uu' (see lLaskarll9 90'). 
It should be noted that in terms of both heliocentric coor- 
dinates or Jacobi coordinates, the indirect part does not 
contribute to the secular system. We have provided here 
the expression of the indirect part for the computation of 
non-secular inequalities. 

6. Conclusion 

We have presented a self-contained exposition of the ex- 
pansion of the three-body Hamiltonian in canonical helio- 
centric coordinates in term of Legendre polynomials at any 
order in the ratio of semi-major axes a, which can also be 
adapted in the case of Jacobi coordinates, where only the 
indirect part differs. We have included here all the neces- 
sary material that allows one to write explicitly the secular 
Hamiltonian at order a 10 for the planar case, a 5 for the spa- 
tial case expressed in terms of the mutual inclination, and 
a 3 for the spatial case in a fixed reference frame. With the 
additional computation of the required Hansen coefficients, 
the expressions of the Tisserand functions F n can also be 
used for a straightforward computation of the expression of 
a non-secular inequality kM + k'M'. 

As the algorithms that are presented here are very sim- 
ple, we have not added to this paper any tables in electronic 
form. Indeed, the reader who needs to use expressions of 



high order that do not appear in the paper, will have no 
problems in exploiting the algorithms given here to derive 
the required expressions. The written results of the paper 
can then be used to check his computations for the lowest 
orders. For example, the computation of -Fiocb m the spatial 
case takes less than 8 minutes on an average 8 core desk- 
top computer in exact rational arithmetics using TRIP, but 
has 2343926 terms, while F50 needs only 7.5 seconds with 
already 164151 terms. As the algorithms require only a few 
lines of code (fewer than 20 in our case), one understands 
that it is preferable to compute the terms when they are 
needed than to store the values in electronic form. 
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Appendix A: Computation of the Hansen coefficients 

The Hansen coefficients X^' m (e) are defined as the Fourier coefficients 



(:)"■ 



J2 X%> m (e)e lkM . (A.l) 



fe=-c 



For arbitrary values of k G Z, the Hansen coefficients X k ' (e) can be expressed in an explicit ma nner as an infinite series 
involving Bessel functions and hypergeometric functions (Hansen 1855; Hill 1875; Tisserand 1899), or generalized Laplace 
coefficients (|Laskarll2005l ). but for k = 0, X£' m (e) reduces to a simple polynomial in e, 1/e, Vl — e 2 , and 1/Vl — e 2 . 
The literature on the computation of the Hansen coefficients has been huge since the original work of lHanse n| p55h . 
but we do not review it here. We concentrate on the obtention of the coefficients Xj , ' m (e) and X^~ n,m (e), for n > and 
< m < n that are required to computate the averaged planetary or lunar Hamiltonian. Using 



r 2 r 



dM = dv = -dE (A.2) 



a 2 v/l 



e 



2 a 



where v is the true anomaly and E the eccentric anomaly, we have for n > 2 







i i r 2 w^ n - 2 



y/l _ e 2 2tt 7 V r 
1 1 /" 27r 

(1 _ e2) „-3/2^y o 



(l + ecosi;)"- 2 e lml, d W (A.3) 



1 [(»-g')/2] (n _ 2) , ^ " 



(1 - e 2 )"- 3 / 2 ^ Z! (m + 1)1 (n-2-(m + 2/))! V 2 

We still need to consider the case n = 1 for which the expansion in true anomaly v is not suitable. We have immediately, 
using (IA.6I) 

A" 1 ' = 1 ; A" 1 ' 1 = V/T 7^' 1 . (A.4) 

It is important to note that from (|A.3|) . we have 

x -n,m = Q for n > 2 and m > n _ ! _ (A.5) 



The Hansen coefficients A n,m are given in Table IA~T1 for n > 2 and m < n — 2. On the other hand, the computation 
of Xq ,m (e) for n > is not as straightforward, as a direct expansion of 



X ° ~27 



1 f 27T fr^ ™ +1 



i ( A -6) 
— / (1 - ecos£:) n+1 " m (cos£: - e + iy/l- e 2 sin E) m dE 
27r Jo 



in eccentric anomaly is far more complicated than the previous expression (IA.3|) . One can still perform some formal 
expansion, but this would not be very useful, as it would contain expressions with many summations that cannot be 
easily reduced to a single summation as in (|A.3[) . In fa ct, the only explicit computation of A"' m (e) available in the 
literature, was obtained through a complex process in (|Hansenlll855t iHiil 118751 : lTisserandlll899[ ). using the auxiliary 
variable ft = (1 — \fl — e 2 )/e. In this case, it can be shown that Xg' m (e) is expressed as a finite hypergeometr ic series in 
(3 2 . U sing a relation among hypergeometric series due to Gauss, and the change of variable e = 2/3/(1 + B 2 ). iTisserandl 
(|1899[) provides a finite expression of Xq ' m (e) in term of hypergeometric series of e 2 . We present here a more direct 
demonstration for the obtention of an explicit expression of X^ m (e) which is a recurrence using the relation among 
Hansen coefficients 

o = n , 9 _ — -(w-ljXo 1 +(n + m)X ' . (A.7) 

This relation appears in (|Hugheslll98TI), but it was shown by iLaskarl (|2005f) that it is equivalent to a recurrence relation 
obtained on Laplace coefficients by lLaplacel (|1785f ). The computation of Xq'° for n > using (|A.6|) is straightforward 
and gives 

[(«+i)/2] / . ,s, / \ 2; 
g SB^WG) ■ < A - 8 » 
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n 


m 




n 


m 


A o ' 


2 





1 


8 





+1 + 15/2e 2 +45/8e 4 + 5/16e 6 


3 





1 


8 


1 


+3e + 15/2e 3 + 15/8e 5 


3 


1 


+l/2e 


8 


2 


+ 15/4 e 2 + 15/4 e 4 + 15/64 e 6 


4 





+ 1 + 1/2 e 2 


8 


3 


+5/2e 3 + 15/16e 5 


4 


1 


+ e 


8 


4 


+ 15/16 e 4 + 3/32 e 6 


4 


2 


+ l/4e 2 


8 


5 


+3/16 e 5 


5 





+ 1 + 3/2e 2 


8 


6 


+1/64 e 6 


5 


1 


+3/2e + 3/8e 3 


9 





+ 1 + 21/2 e 2 + 105/8 e 4 + 35/16 e 6 


5 


2 


+3/4 e 2 


9 


1 


+7/2e + 105/8e 3 + 105/16 e 5 + 35/128 e 7 


5 


3 


+ l/8e 3 


9 


2 


+21/4 e 2 +35/4 e 4 + 105/64 e 6 


6 





+ 1 + 3e 2 + 3/8 e 4 


9 


3 


+35/8e 3 + 105/32e 5 + 21/128e 7 


6 


1 


+2e + 3/2 e 3 


9 


4 


+35/16 e 4 + 21/32 e 6 


6 


2 


+3/2 e 2 + 1/4 e 4 


9 


5 


+21/32 e 5 + 7/128 e 7 


6 


3 


+ 1/2 e 3 


9 


6 


+7/64 e 6 


6 


4 


+ l/16e 4 


9 


7 


+ l/128e 7 


7 





+ 1 + 5e 2 + 15/8 e 4 


10 





+ 1 + 14 e 2 + 105/4 c 4 + 35/4 e 6 + 35/128 e 8 


7 


1 


+5/2e + 15/4e 3 + 5/16e 5 


10 


1 


+4 e + 21e 3 + 35/2 e 5 + 35/16 e 7 


7 


2 


+5/2 e 2 + 5/4 e 4 


10 


2 


+7e 2 + 35/2 e 4 + 105/16 e 6 + 7/32 e 8 


7 


3 


+5/4 e 3 + 5/32 e 5 


10 


3 


+7e 3 + 35/4 e 5 + 21/16 e 7 


7 


4 


+5/16 e 4 


10 


4 


+35/8 e 4 + 21/8 e 6 + 7/64 e 8 


7 


5 


+ 1/32 e 5 


10 


5 


+7/4 e 5 + 7/16 e 7 








10 


6 


+7/16 e 6 + 1/32 e 8 








10 


7 


+l/16e 7 








10 


8 


+1/256 e 8 



Table A.l. X n 



is the polynomial part of the Hansen coefficients X Q n ' m . 



Hansen coefficients are X 
have Xq 1 ' = 1 and Xq' 1 ' 1 



= x -"<"7(i 



2\n-3/2 



For n > 2, we have also X 



For (0 < n < 10,0 < m < n - 2), the full 

n,n — 1 



- A o 



= 0. For n = 1, we 



= (Vl^2 



l)/e. 



The computation of ^ can also be performed using (|A.6[) . Changing E to — E, it is immediate to see that the part in 
sin E cancel, and we are left with 



[n/2] 



X^ = -(n + 2) V 



^ n(l + l)!(n-2l)! V2 



2;+i 



(A.9) 



We can now prove by recurrence using the relation (|A.7|) that the general form of Xq'" 1 for n > and < m < n is 



X^ m = (-iy 



l iii m [(n+l-m)/2] 
(n + 1 + mj ! » ^ 



(n + 1 — m) ! 



(n + 1)! ^ Z!(m + Z)!(n + l-m-2Z)! \2 



m+2/ 



(A.10) 



The computation is delicate but straightforward. One can first show that the two first elements (for / = 0) of the 
polynomial expressions of Xq • m ~ 1 and X^ m ^ 2 in (|A.7|) cancel. We then change the index from I to I + 1 and show that 

the general term of the sum in (|A.7|) gives the general term of Xq'™ 1 , a nd that the last term of Xq'™ 1-2 will give the last 
term of Xq'" 1 . The expression (IA.10[) is equivalent to the expression of (lTisserandl lT899). The expressions of the Hansen 
coefficients Xq'" 1 are given in Table lA~2l for < n < 10 and < m < n. 
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n 


m 




n 


m 


\rn,m 








1 


7 





+ 1 + 14e 2 + 105/4 e 4 + 35/4 e 6 + 35/128 e 8 


1 





+ 1 + 1/2 e 2 


7 


1 


-9/2 e - 189/8 e 3 - 315/16 e 5 - 315/128 e 7 


1 


1 


-3/2 e 


7 


2 


+45/4 e 2 + 225/8 e 4 + 675/64 e 6 +45/128e 8 


2 





+ 1 + 3/2e 2 


7 


3 


-165/8 e 3 - 825/32 e 5 - 495/128 e 7 


2 


1 


-2e - 1/2 e 3 


7 


4 


+495/16 e 4 + 297/16 e 6 + 99/128 e 8 


2 


2 


+5/2 e 2 


7 


5 


-1287/32 e 5 - 1287/128 e 7 


3 





+ l + 3e 2 + 3/8e 4 


7 


6 


+3003/64 e 6 +429/128 e 8 


3 


1 


-5/2 e - 15/8 e 3 


7 


7 


-6435/128 e 7 


3 


2 


+ 15/4e 2 + 5/8 e 4 


8 





+ 1 + 18e 2 + 189/4 e 4 + 105/4 e 6 + 315/128 e 8 


3 


3 


-35/8 e 3 


8 


1 


-5e - 35e 3 - 175/4e 5 - 175/16e 7 - 35/128e 9 


4 





+ 1 + 5e 2 + 15/8 e 4 


8 


2 


+55/4 e 2 + 385/8 e 4 + 1925/64 e 6 + 385/128 e 8 


4 


1 


-3e - 9/2 e 3 - 3/8 e 5 


8 


3 


-55/2 e 3 - 825/16 e 5 - 495/32 e 7 - 55/128 e 9 


4 


2 


+21/4e 2 + 21/8e 4 


8 


4 


+715/16 e 4 + 715/16 e 6 + 715/128 e 8 


4 


3 


-7e 3 - 7/8 e 5 


8 


5 


-1001/16e 5 - 1001/32e 7 - 143/128e 9 


4 


4 


+63/8 e 4 


8 


6 


+5005/64 e 6 + 2145/128 e 8 


5 





+ 1 + 15/2 e 2 + 45/8 e 4 + 5/16 e 6 


8 


7 


-715/8 e 7 - 715/128 e 9 


5 


1 


-7/2 e - 35/4 e 3 - 35/16 e 5 


8 


8 


+ 12155/128 e 8 


5 


2 


+7e 2 + 7e 4 + 7/16 e 6 


9 





+ 1 + 45/2 e 2 + 315/4 e 4 + 525/8 e 6 + 1575/128 e 8 + 63/256 e 10 


5 


3 


-21/2 e 3 - 63/16 e 5 


9 


1 


-11/2 e - 99/2 e 3 - 693/8 e 5 - 1155/32 e 7 - 693/256 e 9 


5 


4 


+ 105/8 e 4 + 21/16 e 6 


9 


2 


+33/2 e 2 + 77e 4 + 1155/16 e 6 + 231/16 e 8 + 77/256 e 10 


5 


5 


-231/16 e 5 


9 


3 


- 143 /4 e 3 - 3003/32 e 5 - 3003/64 e 7 - 1001 /256 e 9 


6 





+ 1 + 21/2 e 2 + 105/8 e 4 + 35/16 e 6 


9 


4 


+ 1001/16e 4 + 3003/32e 6 + 3003/128e 8 + 143/256e 10 


6 


1 


-4e - 15e 3 - 15/2 e 5 - 5/16 e 7 


9 


5 


-3003 /32 e 5 - 5005 /64 e 7 - 2145 /256 e 9 


6 


2 


+9e 2 + 15e 4 +45/16 e 6 


9 


6 


+ 1001/8 e 6 + 429/8 e 8 +429/256 e 10 


6 


3 


-15e 3 - 45/4 e 5 - 9/16 e 7 


9 


7 


-2431/16 e 7 - 7293/256 e 9 


6 


4 


+165/8 e 4 + 99/16 e 6 


9 


8 


+21879/128 e 8 + 2431/256 e 10 


6 


5 


-99/4 e 5 - 33/16 e 7 


9 


9 


-46189/256 e 9 


6 


6 


+429/16 e 6 


10 





+ 1 + 55/2 e 2 + 495/4 e 4 + 1155/8 e 6 + 5775/128 e 8 + 693/256 e 10 








10 


1 


-6e - 135/2e 3 - 315/2e 5 - 1575/16e 7 - 945/64e 9 - 63/256e u 








10 


2 


+39/2 e 2 + 117e 4 + 2457/16 e 6 + 819/16 e 8 + 819/256 e 10 








10 


3 


-91/2e 3 - 637/4e 5 - 1911/16e 7 - 637/32e 9 - 91/256e n 








10 


4 


+1365/16 e 4 + 5733/32 e 6 + 9555/128 e 8 + 1365/256 e 10 








10 


5 


-273/2 e 5 - 1365/8 e 7 - 585/16 e 9 - 195/256 e 11 








10 


6 


+ 1547/8e 6 + 1105/8e 8 +3315/256e 10 








10 


7 


-1989/8 e 7 - 5967/64 e 9 - 663/256 6" 








10 


8 


+37791/128 e 8 + 12597/256 e 10 








10 


9 


-20995/64 e 9 - 4199/256 e 11 








10 


10 


+88179/256 e 10 



Table A. 2. Hansen coefficients Xq'™ 1 for (0 < n < 10, < m < n). 



Appendix B: Computation of Tisserand polynomials and Tisserand functions 

Expansions valid for all inclinations were already given by lHanse n (1855), and iTisserandl f|1885h in his researches on 
asteroidal motions, but these derivations are rather complex. Although in (|B.3|) we present some of Tisserand's compu- 
tations, we make first, as in the planar case some direct expansions and reordering of terms to obtain expressions that 
are as compact as possible. 

B.l. Tisserand functions with respect to mutual inclination 

With x = u — u', y = u + u',fi = cos 2 ( J/2), v = sin 2 ( J/2), we have 

uu' = /.icosx + ^cos y , (B.l) 
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and with a straightforward expansion in complex notation, we obtain 
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which we can reorder as 



with 



ml n ki+k2 7 k 3 +k 4 

K 1 + « 2 + K 3 + K 4 — "I 

0<fc; <m 



(u-u') m - ^rE E CiS^e-V^-ri (B.3) 

r— p=—r 



rain(r,m— r-p) ,,2Z+p„m-2Z-p 

C( m ) = V - - 



HQ + p)\(r - l)\(m - p - r - l)\ ' 

i=max(0,-p) V 1 J K ' K y ' 



With the expansion (|2.8p . we have then 



with 



[n/2] n-1kn~1k-r 

F„=E^E E C^^e-'e* (B.5) 

fc=0 r=0 p= — r 



[n-2k)\ _ {-If (2n-2fc)! 



We can then reorder (IB. 51) with s = fc + r as 



n min(s,n — s) n — fc — s 

F »=E E E ffn.fcC^-^e-^^- 3 '-*') , (B.7) 

s— fc— p—k — s 



exchange the sum in p and fc 



n n-s min(s+p,n-s—p,s,n—s) 

fn=EE E C_ 2 1 P e'^e^ C-»-»0 (B.8) 

s=0 p=— s fc=0 



and set q = p + s which gives 



n n I min(q,n—q,s,n—s) 

F " = EE E WW,-. I e^C'-Je 1 * . (B.9) 



s=0 q=0 \ k=0 

That is, with x — u — u' and y = u + u' , 



with 



that is 



F n = EE 2 ^ ) ^' i/ ) el " ( "~ 2s)eW ' (n ~ 29) ( R1 °) 

s=0 q=Q 



min(g,n— q,s,n — s) 

Q^(J)= E ?»,*<fc,nV.> (B.ll) 

fc=0 

mln(q,n-q,s,n-s) min(s-fc,n-fc- 9 ) 2i+ 9 -s n-2k-q+s-2l 

efc>(-o = E E i!(f+g _ / a)!(j _ A; _ , (ra _ fc _ g _ , ■ (B.12) 

fc=0 I=max(0,s- 5 ) V ^ ^ V ^ V. 1 ) 
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B.2. Tisserand functions with respect to a fixed reference frame 

In the case of a fixed reference frame fSect !4.2j) . the expressions differ slightly, but are more complicated. We have (|4.13l) 

uu' = acosx — a' smx + bcosy — b' s'my , (B.13) 

that is, with /i* = a + iaf, v* = b + ib', 

u-u' = i [^e lx + Jl,e- lx + + V,e~ iy ] . (B.14) 



This is similar to the previous expression (|B.1[) , and the computation of the Tisserand functions F n will follow the same 
lines and gives the more general expression 

n n 

F n =Y J J2 2«?<r e ™ ( ™~ 2s) e m ' {n ~ 2q) (B.15) 

s=0 q=0 

with 

min(q,n-q,s,n-s) min(s-fe,n-fc-g) l+q-s — I n-k-q-l — s-k-l 

a£Sfe,«.)- E E l!(l+ ^ e) C-l - t -,-or <B - 16) 

fe=0 i=max(Q,s-g) V y ; V ; V * ' 

The coefficients Q^g (J) (|B.12[> and Qs™q (/U*, i/*) (|B.16|I in the Tisserand functions are very similar. Indeed, using /i = 1 — z/ 
and |/x*| = 1 — (|4.12l) . we have 

q^o = M'-^-'-'^-U-fl-M. s&V*^*) - Arx- 9 -*4-.,«- a -.(M)> (B.17) 

with z; = 1 1/, | = sin 2 (J/2) and 

roin(«,n-g,s,n-s) min(s-fc,n-fc-ij) _ _ x \2l x —2k+2s—2l 

= E ; _ ?) + g _,)!(/_ fe - !(n-ife-,-0! ' ^ 

Tisserand simplification 

The above expression (|B.18[) is a double sum. In the mutual inclination case, iTisseraiicfl Jl885) could reduce it to a single 
sum using hypergeometric functions. He first considers the differential equation satisfied by Legendre polynomials 

(l-z 2 )^i _ 2 ^+n(n + l)P n =0 . (B.19) 

As i 7 ^ = P n {z) CEa l2.8[) with z = cos a; + z/(cosy — cosir) (Eg. IB. II) . we have from (IB. 19ft . 

B 2 F BF 1 B 2 F ~\ B 2 F 

1/(1 - i/)^ + (1 - 2z/& + . + - + n(n + l)F n = . (B.20) 

Bv A Bv 1 — v Bx z v By 1 

Using expressions (jB.101 [B~TT|) . one then replaces F n in (1B.20[) by 

F » = E E(! - "^^(^ , (B.21) 
fc i 

which leads to 

££(1 - u)V^ Me ^» = , (B.22) 
fe z 

where 



4"'M = 1/(1 - v )—^~ + t 1 + 2/ - 2 ( fc + 1 + !H + (n - fc - l)(n + k + l + 1)A%{ . (B.23) 

The equality (|B.22|) must be satisfied for all x and y, thus all the coefficients A^}{v) (|B.23j) are equal to 0. The solutions 
of |B33l) are 

A^J{v) = K ( ^F(k + l-n,k + l + n+l;l + 2l;u) (B.24) 

or equivalently 

= K^\ n _ g _ s F(-2s, 2n -2s + l;2n-2q-2s+ 1; i/) , (B.25) 
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where Kuj is an unknown constant and F an hypergeometric function (e.g. IWhittaker fc Watsonlll927t ). Let us assume 
that the quantity 2n—2q—2s+l is positive. If it is not the case, one can make the change of variable (s, q) — > [n—s 1 , n—q'). 
Then 2n - 2q' - 2s' + 1 is positive and since F n is real, Qs$(i>) = 2^/0) (|B7T0]) . From fR25]l . one thus has 

,(„) M = ^(«) (2 a )!(2n-2g-2 S )! ^ fc (2n-2s + fc)! 

9 -.,n-9-»W 9-»,n-9-. (2 n - 2s)! ^ ' (2s - fr)!(2n - 2q - 2s + fc)! fc! ' 1 ' j 

iTisseraiicfl ()1885() needs then some lengthy computations to determine K^ 1 ) from the coefficient of v n in F n . Here, we use 
the expression (|B~T2l) with (|BTT7| and ([H726j) . With fj, = 1 - i/, we get 

(2n-2 g -2s)!(2n)! (2n)! 1 

a -.,n- g -.V ; (2n-2s)!(2n-2g)! V ' 2 2n n\ + q - s)\(s - l)l(n - q - 1)1 ' 1 ' ; 



Calculating the coefficient of the term x s in (1 + x) n q (l + x) q = (1 + x) n , one finds 

1 _ n\ 

l\(l + q-s)\{s-l)\{n-q-l)\ ~ s!(n - s)\q\{n - q)\ 



min(s,n — q) ^ ^ 

51 Wl 4-n- - - n -IV ~ S \(r, - S )\n\(r, - „V ' (B.28) 

i=max(0,s— q) 



Thus, from (|B.27[) , we obtain 



K (n) = 1 (2n-2s)!(2n-2 g )! 

q -s,n- q -s 2 2n s\{n - s)\q\{n - q)\{2n - 2q - 2s)\ ' 1 ' ' 



Finally, the most general case, where inclinations are defined with respect to a fixed reference plane, can be as well 
derived from (|B~T71) and (|R29|t . We have 

n n 
s=0 q=0 



where 



and 



Q$Q**,v.) = { (B.31) 
^r q K s+q - n A ( ^ s+q -M) else, 



(n) 1 (2s)!(2n-2g)! A fc (2n-2s + fc)! x fc 

,-.,„-,-. W 2 2 ™s!(n-s)!g!(n-g)! ^ l J (2s - fc)!(2n - 2q - 2s + fc)! fc! 

In (|B.31[) . we used the fact that is the complex conjugate of Q^!: s „_„ since F„ is real. 



